
####################################################
####          Replication Code for              ####
####     Spatial Interdependence                ####
####     and Instrumental Variable Models       ####
####################################################
####                                            ####
####                9/ 06/2018                  ####
####  this file runs rcode to create            ####
####            Figures in Appendix             ####
####################################################




#### First Figures for Appendix A, i.e., with misspecification in W
#### load simulation output
load("simulation_output_RR.rda")
summary(dat)



### some housekeeping: we change variable names
dat$Method = revalue(dat$Method, c("2sls-robust" = "2slsRobust"))
unique(dat$Method)

### We subset the output to drop 2sls with robust standard errors as these are not analyzed
### No substantial improvement with robust SE
dat = dat[dat$Method %in% c("ols","2sls","SAR-iv"),]
### Next we subset the output to include only the simulations with correlation equal to 0.5, i.e. high level of non-spatial endogeneity
dat = dat[dat$cor == 0.5, ]

### calculate coverage statistic, absolut and squared error
dat$coverage = 0
dat$coverage = ifelse(dat$treat <= (dat$Coef + (qnorm(0.975)*dat$SE)) & dat$treat >= (dat$Coef - (qnorm(0.975)*dat$SE)), 1, dat$coverage)
dat$abs.error = abs(dat$treat - dat$Coef)
dat$sqr.error = (dat$treat - dat$Coef)^2



### for some simulated datasets with misspecification in W, the spatial model did not estimate. We drop those simulation runs for all estimation methods, it does not make a difference when we do not drop the corresponding estimations for LM and standard IV methods
missingSpatial <- which(is.na(dat$Coef)==T)
### delete corresponding LM and IV
dat2 <- dat[-c((missingSpatial - 1), (missingSpatial - 2)), ]





##### subset output for strong IV and large N setting
data_iv_strong = subset(dat, iv.coef == 1.5 & N == 200)
#### data summary
data_plot = ddply(data_iv_strong, .(Method, rho.y, rho.z, misspec), summarize, meanAbs = mean(abs.error, na.rm = T), medianAbs = median(abs.error, na.rm = T), coverage = mean(coverage, na.rm = T),  rmse = sqrt(mean(sqr.error, na.rm = T)))


### Rename spatial 2sls for the plot
### Create Variable that will serve as Label for the Row
data_plot$Method = revalue(data_plot$Method, c("SAR-iv" = "s-2sls"))
data_plot$misspecification = case_when(data_plot$misspec == 0 ~ "correct \n W", data_plot$misspec == 0.5 ~ "misspecification \n probability = 0.5", data_plot$misspec == 1 ~"wrong \n W")




p = ggplot(data = data_plot, aes(y = medianAbs, x = rho.y))
p = p + geom_point(aes(colour = Method, shape = Method), size = 4.5, alpha = 1, position=position_dodge(width = 0.1))
p = p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p = p + labs(title = "", y="Median Absolute Error", x = expression(rho[y]))+ theme(axis.text.x=element_text(size=13, face = "bold")) + theme(axis.text.y=element_text(size=15, face = "bold")) + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20))
p = p + scale_colour_brewer(palette = "Set1",name = "Estimation Method") + scale_shape_manual(name = "Estimation Method", values = c(15, 17, 16))
p = p + theme(legend.position = "none",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 15), legend.title = element_text(size = 15), panel.border = element_rect(colour = "black", fill=NA, size=2.5))
p = p + facet_grid(misspecification ~ rho.z, labeller = label_bquote(cols = rho[z] ~ "=" ~ .(rho.z))) +  theme(strip.background = element_blank(), strip.placement = "outside") + theme(strip.text.x = element_text(size=14, face = "bold"), strip.text.y = element_text(size=12, face="bold", angle = 45)) + scale_x_continuous(limits = c(-0.1,0.7), breaks = c(0,0.3,0.6),labels = c(0,0.3,0.6))
p = p + theme(legend.position = "bottom", legend.direction = "horizontal",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 24), legend.title = element_text(size = 24), panel.border = element_rect(colour = "black", fill=NA, size=3.5))
p = p + theme(panel.spacing = unit(0.95, "lines")) +  scale_y_continuous(limits = c(0,1.25), breaks = c(0,0.25, 0.5, 0.75, 1, 1.25))
plot(p)
ggsave("FigureA1.pdf", width = 6 * 1.618, height = 6)


##same plot for coverage
p = ggplot(data = data_plot, aes(y = coverage, x = rho.y))
p = p + geom_point(aes(colour = Method, shape = Method), size = 4.5, alpha = 1, position=position_dodge(width = 0.1))
p = p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p = p + labs(title = "", y="Coverage", x = expression(rho[y]))+ theme(axis.text.x=element_text(size=13, face = "bold")) + theme(axis.text.y=element_text(size=15, face = "bold")) + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20))
p = p + scale_colour_brewer(palette = "Set1",name = "Estimation Method") + scale_shape_manual(name = "Estimation Method", values = c(15, 17, 16))
p = p + theme(legend.position = "none",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 15), legend.title = element_text(size = 15), panel.border = element_rect(colour = "black", fill=NA, size=2.5))
p = p + facet_grid(misspecification~ rho.z, labeller = label_bquote(cols = rho[z] ~ "=" ~ .(rho.z))) +  theme(strip.background = element_blank(), strip.placement = "outside") + theme(strip.text.x = element_text(size=14, face = "bold"), strip.text.y = element_text(size=12, face="bold", angle = 45)) + scale_x_continuous(limits = c(-0.1,0.7), breaks = c(0,0.3,0.6),labels = c(0,0.3,0.6))
p = p + theme(legend.position = "bottom", legend.direction = "horizontal",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 24), legend.title = element_text(size = 24), panel.border = element_rect(colour = "black", fill=NA, size=3.5)) +  scale_y_continuous(limits = c(0,1.25), breaks = c(0,0.5,0.95))
p = p + theme(panel.spacing = unit(0.95, "lines")) + geom_hline(aes(yintercept = 0.95),col = "red")
plot(p)
ggsave("FigureA2.pdf", width = 6 * 1.618, height = 6)




#### Figures for Appendix B and C
#### reload data to analyze sims without misspecification for Appendix B & C
#### load simulation output
rm(dat)
load("simulation_output_RR.rda")
summary(dat)



### some housekeeping: we change variable names
dat$Method = revalue(dat$Method, c("2sls-robust" = "2slsRobust"))
unique(dat$Method)

### We subset the output to drop 2sls with robust standard errors as these are not analyzed
### No subtantial improvement with robust SE
dat = dat[dat$Method %in% c("ols","2sls","SAR-iv"),]
### Next we subset the output to include only the simulations without misspecification of W
dat = dat[dat$misspec == 0, ]

### calculate coverage statistic, absolut and squared error
dat$coverage = 0
dat$coverage = ifelse(dat$treat <= (dat$Coef + (qnorm(0.975)*dat$SE)) & dat$treat >= (dat$Coef - (qnorm(0.975)*dat$SE)), 1, dat$coverage)
dat$abs.error = abs(dat$treat - dat$Coef)
dat$sqr.error = (dat$treat - dat$Coef)^2




########################################################################################################################### ###########################################################################################################################
###########################################################################################################################
###########################################################################################################################
#### Now we create the same plots under different scenarios
#### Weaker IV, N = 200
#### Figure B3 and B4

### Subset data
data_iv_med = subset(dat, iv.coef == 0.75 & N == 200)

#### Calculate summary stats for Figure
data_plot = ddply(data_iv_med, .(Method, rho.y, rho.z, cor), summarize, meanAbs = mean(abs.error), medianAbs = median(abs.error), coverage = mean(coverage),  mse = mean(sqr.error))

### Rename spatial 2sls for the plot
### Create Variable that will serve as Label for the Row
data_plot$Method = revalue(data_plot$Method, c("SAR-iv" = "s-2sls"))
data_plot$nonspat = case_when(data_plot$cor == -0.5 ~ "non-spatial \n endogeneity = -0.5", data_plot$cor == 0 ~ "non-spatial \n endogeneity = 0", data_plot$cor == 0.5 ~"non-spatial \n endogeneity = 0.5")


p = ggplot(data = data_plot, aes(y = medianAbs, x = rho.y))
p = p + geom_point(aes(colour = Method, shape = Method), size = 4.5, alpha = 1, position=position_dodge(width = 0.1))
p = p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p = p + labs(title = "", y="Median Absolute Error", x = expression(rho[y]))+ theme(axis.text.x=element_text(size=13, face = "bold")) + theme(axis.text.y=element_text(size=15, face = "bold")) + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20))
p = p + scale_colour_brewer(palette = "Set1",name = "Estimation Method") + scale_shape_manual(name = "Estimation Method", values = c(15, 17, 16))
p = p + theme(legend.position = "none",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 15), legend.title = element_text(size = 15), panel.border = element_rect(colour = "black", fill=NA, size=2.5))
p = p + facet_grid(nonspat~ rho.z, labeller = label_bquote(cols = rho[z] ~ "=" ~ .(rho.z))) +  theme(strip.background = element_blank(), strip.placement = "outside") + theme(strip.text.x = element_text(size=14, face = "bold"), strip.text.y = element_text(size=12, face="bold", angle = 45)) + scale_x_continuous(limits = c(-0.1,0.7), breaks = c(0,0.3,0.6),labels = c(0,0.3,0.6))
p = p + theme(legend.position = "bottom", legend.direction = "horizontal",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 24), legend.title = element_text(size = 24), panel.border = element_rect(colour = "black", fill=NA, size=3.5))
p = p + theme(panel.spacing = unit(0.95, "lines")) +  scale_y_continuous(limits = c(0,1.25), breaks = c(0,0.25, 0.5, 0.75, 1, 1.25))
plot(p)
ggsave("FigureB3.pdf", width = 6 * 1.618, height = 6)


##same plot for coverage
p = ggplot(data = data_plot, aes(y = coverage, x = rho.y))
p = p + geom_point(aes(colour = Method, shape = Method), size = 4.5, alpha = 1, position=position_dodge(width = 0.1))
p = p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p = p + labs(title = "", y="Coverage", x = expression(rho[y]))+ theme(axis.text.x=element_text(size=13, face = "bold")) + theme(axis.text.y=element_text(size=15, face = "bold")) + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20))
p = p + scale_colour_brewer(palette = "Set1",name = "Estimation Method") + scale_shape_manual(name = "Estimation Method", values = c(15, 17, 16))
p = p + theme(legend.position = "none",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 15), legend.title = element_text(size = 15), panel.border = element_rect(colour = "black", fill=NA, size=2.5))
p = p + facet_grid(nonspat~ rho.z, labeller = label_bquote(cols = rho[z] ~ "=" ~ .(rho.z))) +  theme(strip.background = element_blank(), strip.placement = "outside") + theme(strip.text.x = element_text(size=14, face = "bold"), strip.text.y = element_text(size=12, face="bold", angle = 45)) + scale_x_continuous(limits = c(-0.1,0.7), breaks = c(0,0.3,0.6),labels = c(0,0.3,0.6))
p = p + theme(legend.position = "bottom", legend.direction = "horizontal",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 24), legend.title = element_text(size = 24), panel.border = element_rect(colour = "black", fill=NA, size=3.5)) +  scale_y_continuous(limits = c(0,1.25), breaks = c(0,0.5,0.95))
p = p + theme(panel.spacing = unit(0.95, "lines")) + geom_hline(aes(yintercept = 0.95),col = "red")
plot(p)
ggsave("FigureB4.pdf", width = 6 * 1.618, height = 6)
##same plot for coverage





########################################################################################################################### ###########################################################################################################################
###########################################################################################################################
###########################################################################################################################
#### Now we create the same plots under different scenarios
#### strong IV, N = 50
#### Figure B5 and B6


data_iv = subset(dat, iv.coef == 1.5 & N == 50)
data_plot = ddply(data_iv, .(Method, rho.y, rho.z, cor ), summarize, meanAbs = mean(abs.error), medianAbs = median(abs.error), coverage = mean(coverage),  mse = mean(sqr.error))


### Rename spatial 2sls for the plot
### Create Variable that will serve as Label for the Row
data_plot$Method = revalue(data_plot$Method, c("SAR-iv" = "s-2sls"))
data_plot$nonspat = case_when(data_plot$cor == -0.5 ~ "non-spatial \n endogeneity = -0.5", data_plot$cor == 0 ~ "non-spatial \n endogeneity = 0", data_plot$cor == 0.5 ~"non-spatial \n endogeneity = 0.5")

p = ggplot(data = data_plot, aes(y = medianAbs, x = rho.y))
p = p + geom_point(aes(colour = Method, shape = Method), size = 4.5, alpha = 1, position=position_dodge(width = 0.1))
p = p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p = p + labs(title = "", y="Median Absolute Error", x = expression(rho[y]))+ theme(axis.text.x=element_text(size=13, face = "bold")) + theme(axis.text.y=element_text(size=15, face = "bold")) + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20))
p = p + scale_colour_brewer(palette = "Set1",name = "Estimation Method") + scale_shape_manual(name = "Estimation Method", values = c(15, 17, 16))
p = p + theme(legend.position = "none",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 15), legend.title = element_text(size = 15), panel.border = element_rect(colour = "black", fill=NA, size=2.5))
p = p + facet_grid(nonspat~ rho.z, labeller = label_bquote(cols = rho[z] ~ "=" ~ .(rho.z))) +  theme(strip.background = element_blank(), strip.placement = "outside") + theme(strip.text.x = element_text(size=14, face = "bold"), strip.text.y = element_text(size=12, face="bold", angle = 45)) + scale_x_continuous(limits = c(-0.1,0.7), breaks = c(0,0.3,0.6),labels = c(0,0.3,0.6))
p = p + theme(legend.position = "bottom", legend.direction = "horizontal",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 24), legend.title = element_text(size = 24), panel.border = element_rect(colour = "black", fill=NA, size=3.5)) +  scale_y_continuous(limits = c(0,1), breaks = c(0,0.25, 0.5, 0.75, 1))
p = p + theme(panel.spacing = unit(0.95, "lines"))
plot(p)
ggsave("FigureB5.pdf", width = 6 * 1.618, height = 6)


##same plot for coverage
p = ggplot(data = data_plot, aes(y = coverage, x = rho.y))
p = p + geom_point(aes(colour = Method, shape = Method), size = 4.5, alpha = 1, position=position_dodge(width = 0.1))
p = p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p = p + labs(title = "", y="Coverage", x = expression(rho[y]))+ theme(axis.text.x=element_text(size=13, face = "bold")) + theme(axis.text.y=element_text(size=15, face = "bold")) + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20))
p = p + scale_colour_brewer(palette = "Set1",name = "Estimation Method") + scale_shape_manual(name = "Estimation Method", values = c(15, 17, 16))
p = p + theme(legend.position = "none",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 15), legend.title = element_text(size = 15), panel.border = element_rect(colour = "black", fill=NA, size=2.5))
p = p + facet_grid(nonspat~ rho.z, labeller = label_bquote(cols = rho[z] ~ "=" ~ .(rho.z))) +  theme(strip.background = element_blank(), strip.placement = "outside") + theme(strip.text.x = element_text(size=14, face = "bold"), strip.text.y = element_text(size=12, face="bold", angle = 45)) + scale_x_continuous(limits = c(-0.1,0.7), breaks = c(0,0.3,0.6),labels = c(0,0.3,0.6))
p = p + theme(legend.position = "bottom", legend.direction = "horizontal",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 24), legend.title = element_text(size = 24), panel.border = element_rect(colour = "black", fill=NA, size=3.5)) +  scale_y_continuous(limits = c(0,1.25), breaks = c(0,0.5,0.95))
p = p + theme(panel.spacing = unit(0.95, "lines")) + geom_hline(aes(yintercept = 0.95),col = "red")
plot(p)
ggsave("FigureB6.pdf", width = 6 * 1.618, height = 6)






##### next we create plots for Appendix C,
#### Figure C.7
#### strong IV, large N
data_iv_strong = subset(dat, iv.coef == 1.5 & N == 200)
#### calculate summaries
data_plot = ddply(data_iv_strong, .(Method, rho.y, rho.z, cor ), summarize, meanAbs = mean(abs.error), medianAbs = median(abs.error), coverage = mean(coverage),  rmse = sqrt(mean(sqr.error)))


### Rename spatial 2sls for the plot
### Create Variable that will serve as Label for the Row
data_plot$Method = revalue(data_plot$Method, c("SAR-iv" = "s-2sls"))
data_plot$nonspat = case_when(data_plot$cor == -0.5 ~ "non-spatial \n endogeneity = -0.5", data_plot$cor == 0 ~ "non-spatial \n endogeneity = 0", data_plot$cor == 0.5 ~"non-spatial \n endogeneity = 0.5")


p = ggplot(data = data_plot, aes(y = rmse, x = rho.y))
p = p + geom_point(aes(colour = Method, shape = Method), size = 4.5, alpha = 1, position=position_dodge(width = 0.1))
p = p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p = p + labs(title = "", y="RMSE", x = expression(rho[y]))+ theme(axis.text.x=element_text(size=13, face = "bold")) + theme(axis.text.y=element_text(size=15, face = "bold")) + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20))
p = p + scale_colour_brewer(palette = "Set1",name = "Estimation Method") + scale_shape_manual(name = "Estimation Method", values = c(15, 17, 16))
p = p + theme(legend.position = "none",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 15), legend.title = element_text(size = 15), panel.border = element_rect(colour = "black", fill=NA, size=2.5))
p = p + facet_grid(nonspat~ rho.z, labeller = label_bquote(cols = rho[z] ~ "=" ~ .(rho.z))) +  theme(strip.background = element_blank(), strip.placement = "outside") + theme(strip.text.x = element_text(size=14, face = "bold"), strip.text.y = element_text(size=12, face="bold", angle = 45)) + scale_x_continuous(limits = c(-0.1,0.7), breaks = c(0,0.3,0.6),labels = c(0,0.3,0.6))
p = p + theme(legend.position = "bottom", legend.direction = "horizontal",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 24), legend.title = element_text(size = 24), panel.border = element_rect(colour = "black", fill=NA, size=3.5))
p = p + theme(panel.spacing = unit(0.95, "lines")) +  scale_y_continuous(limits = c(0,1.25), breaks = c(0,0.25, 0.5, 0.75, 1, 1.25))
plot(p)
ggsave("FigureC7.pdf", width = 6 * 1.618, height = 6)

#### Figure C.8
#### weak IV, large N
data_iv = subset(dat, iv.coef == 0.75 & N == 200)
#### calculate summaries
data_plot = ddply(data_iv, .(Method, rho.y, rho.z, cor ), summarize, meanAbs = mean(abs.error), medianAbs = median(abs.error), coverage = mean(coverage),  rmse = sqrt(mean(sqr.error)))


### Rename spatial 2sls for the plot
### Create Variable that will serve as Label for the Row
data_plot$Method = revalue(data_plot$Method, c("SAR-iv" = "s-2sls"))
data_plot$nonspat = case_when(data_plot$cor == -0.5 ~ "non-spatial \n endogeneity = -0.5", data_plot$cor == 0 ~ "non-spatial \n endogeneity = 0", data_plot$cor == 0.5 ~"non-spatial \n endogeneity = 0.5")


p = ggplot(data = data_plot, aes(y = rmse, x = rho.y))
p = p + geom_point(aes(colour = Method, shape = Method), size = 4.5, alpha = 1, position=position_dodge(width = 0.1))
p = p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p = p + labs(title = "", y="RMSE", x = expression(rho[y]))+ theme(axis.text.x=element_text(size=13, face = "bold")) + theme(axis.text.y=element_text(size=15, face = "bold")) + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20))
p = p + scale_colour_brewer(palette = "Set1",name = "Estimation Method") + scale_shape_manual(name = "Estimation Method", values = c(15, 17, 16))
p = p + theme(legend.position = "none",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 15), legend.title = element_text(size = 15), panel.border = element_rect(colour = "black", fill=NA, size=2.5))
p = p + facet_grid(nonspat~ rho.z, labeller = label_bquote(cols = rho[z] ~ "=" ~ .(rho.z))) +  theme(strip.background = element_blank(), strip.placement = "outside") + theme(strip.text.x = element_text(size=14, face = "bold"), strip.text.y = element_text(size=12, face="bold", angle = 45)) + scale_x_continuous(limits = c(-0.1,0.7), breaks = c(0,0.3,0.6),labels = c(0,0.3,0.6))
p = p + theme(legend.position = "bottom", legend.direction = "horizontal",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 24), legend.title = element_text(size = 24), panel.border = element_rect(colour = "black", fill=NA, size=3.5))
p = p + theme(panel.spacing = unit(0.95, "lines")) +  scale_y_continuous(limits = c(0,1.25), breaks = c(0,0.25, 0.5, 0.75, 1, 1.25))
plot(p)
ggsave("FigureC8.pdf", width = 6 * 1.618, height = 6)


#### Figure C.9
#### strong IV, small N
data_iv = subset(dat, iv.coef == 1.5 & N == 50)
#### calculate summaries
data_plot = ddply(data_iv, .(Method, rho.y, rho.z, cor ), summarize, meanAbs = mean(abs.error), medianAbs = median(abs.error), coverage = mean(coverage),  rmse = sqrt(mean(sqr.error)))


### Rename spatial 2sls for the plot
### Create Variable that will serve as Label for the Row
data_plot$Method = revalue(data_plot$Method, c("SAR-iv" = "s-2sls"))
data_plot$nonspat = case_when(data_plot$cor == -0.5 ~ "non-spatial \n endogeneity = -0.5", data_plot$cor == 0 ~ "non-spatial \n endogeneity = 0", data_plot$cor == 0.5 ~"non-spatial \n endogeneity = 0.5")


p = ggplot(data = data_plot, aes(y = rmse, x = rho.y))
p = p + geom_point(aes(colour = Method, shape = Method), size = 4.5, alpha = 1, position=position_dodge(width = 0.1))
p = p + theme_bw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p = p + labs(title = "", y="RMSE", x = expression(rho[y]))+ theme(axis.text.x=element_text(size=13, face = "bold")) + theme(axis.text.y=element_text(size=15, face = "bold")) + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20))
p = p + scale_colour_brewer(palette = "Set1",name = "Estimation Method") + scale_shape_manual(name = "Estimation Method", values = c(15, 17, 16))
p = p + theme(legend.position = "none",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 15), legend.title = element_text(size = 15), panel.border = element_rect(colour = "black", fill=NA, size=2.5))
p = p + facet_grid(nonspat~ rho.z, labeller = label_bquote(cols = rho[z] ~ "=" ~ .(rho.z))) +  theme(strip.background = element_blank(), strip.placement = "outside") + theme(strip.text.x = element_text(size=14, face = "bold"), strip.text.y = element_text(size=12, face="bold", angle = 45)) + scale_x_continuous(limits = c(-0.1,0.7), breaks = c(0,0.3,0.6),labels = c(0,0.3,0.6))
p = p + theme(legend.position = "bottom", legend.direction = "horizontal",legend.key = element_blank(), legend.background = element_rect(colour = "white"), legend.text = element_text(size = 24), legend.title = element_text(size = 24), panel.border = element_rect(colour = "black", fill=NA, size=3.5))
p = p + theme(panel.spacing = unit(0.95, "lines")) +  scale_y_continuous(limits = c(0,1.25), breaks = c(0,0.25, 0.5, 0.75, 1, 1.25))
plot(p)
ggsave("FigureC9.pdf", width = 6 * 1.618, height = 6)
